%% 优化的基本算法理论The basic algorithm theory of optimization
%% 1.基础最优性理论The basic theory of optimization
%%  1.1.最优化问题解的存在性：
% 在实际应用中优化问题的形式多种多样，给定一类具体的优化问题时首先要分析其解
% 的存在性，当优化问题解存在时再考虑设计算法求出问题较优解或最优解。
% Weierstrass定理：设f(x)是定义在[a,b]上连续函数，则存在多项式函数列{fn(x)},
% 使得fn(x)一致收敛于f(x)。由Weierstrass定理可知，定义在有限集上的连续函数一
% 定存在最大或最小值。在实际问题中，定义域可能不是有限集，目标函数可能不连续，
% 因此将Weierstrass定理推广：
% 对于优化问题：min f(x);s.t.xmin<=x<=xmax;,f(x)为适当且闭的函数，假设下面三
% 个条件中任意一个成立：
% (1)f有界;(2)存在常数r，使得下水平集{x|f(x)<=r}非空有界;(3)对于在定义域中任
% 意满足lim(||xk||,'right')=+inf的点列xk，有lim(f(xk))=+inf。
% 则优化问题的解是非空且有限的。
% 当优化问题中f(x)为强拟凸函数时，即函数满足：定义域为凸集，定义域内任意两点
% 之间线段上的函数值小于线段任一端点处函数值时，优化问题存在唯一最优解。
% 对于一般定义在有限凸集上的闭强凸函数，其最优解唯一存在；对于凸函数，其最优
% 解存在。
% matlab中符号函数有关极限，积分，求导的运算如下：
clear,clc
syms x
f1=x^3;a=0;
a1=limit(f1,x,a); %求f关于自变量x在a处的极限值
a2=limit(f1,x,a,'right'); %求f关于自变量x在a处的右极限值
a3=limit(f1,x,a,'left'); %求f关于自变量x在a处的左极限值
f1=exp(-abs(x));
a4=limit(f1,x,-inf); %求f关于自变量x左无穷处极限值
a5=limit(f1,x,inf); %求f关于自变量x右无穷处极限值
f1=x^2+3*x^3;n=2;a=0;b=3;
a6=diff(f1,x,n); %求f关于自变量x的n阶偏导数
a7=int(f1,x,a,b); %求f关于自变量x在区间[a,b]上的定积分
a8=int(f1,x); %求f关于自变量x的不定积分(结果中没有任意常数）
% matlab积分的数值计算方法：
% 匿名函数积分(最大3维)
clear,clc
f1=@(x1)x1.^2;
f2=@(x1,x2)x1.^2+x2.^2;
f3=@(x1,x2,x3)x1.^2+x2.^2+x3.^2;
a=0;b=1;c=0;d=1;e=0;f=1;% 积分上下限
a1=integral(f1,a,b); % 求一元匿名函数f在[a,b]上的积分值,使用
% integral(f,a,b,'RelTol',1e-20)可获得高精度的值
%对于数值向量可以用integral(f,a,b,'RelTol',1e-20,'ArrayValued',true)求积分
a2=integral2(f2,a,b,c,d); %求二元匿名函数f在[a,b]*[c,d]上的二重积分值
a3=integral3(f3,a,b,c,d,e,f); %求三元匿名函数f在[a,b]*[c,d]*[e,f]上的三重积分值
% 数值向量函数积分方法
clear,clc
xmin=0;xmax=10;dx=1e-2;x=xmin:dx:xmax;f=x.^2;% 初始函数
Q=trapz(x,f);% 利用梯形面积微元计算积分
% 利用数值计算方法设计多元数量值函数数值积分：
clear,clc
f=@(x1,x2,x3,x4)x1+x2+x3+x4;
xmin=[0,0,0,0];% 积分下限
xmax=[1,1,1,1];% 积分上限
dx=[1e-1,1e-1,1e-2,1e-2];% 自变量各维度积分微元
dv=prod(dx);% 积分"体积"微元
a=0;
for i1=xmin(1):dx(1):xmax(1)
    for i2=xmin(2):dx(2):xmax(2)
        for i3=xmin(3):dx(3):xmax(3)
            for i4=xmin(4):dx(4):xmax(4)
                a=a+f(i1,i2,i3,i4)*dv;
            end
        end
    end
end
%%  1.2.无约束可微问题的最优性理论：
% 本小节讨论目标函数f(x)连续可微的情况。
% 定义梯度算子gradient(f(x))，则对于可微函数f，若存在向量d满足：
% gradient(f(x))'*d<0，则称d为f在点x处的一个下降方向。
% 若d为f在点x处下降方向，则对于任意T>0，存在0<=t<=T，使得：
% f(x+t*d)<f(x)。
% 一阶最优性条件：若f在定义域可微，x0为局部极小解，则gradient(f(x0))=0。
% 称在定义域中满足gradient(f(x))=0的点x为f的稳定点(驻点或临界点)。
% 二阶最优性条件：设f在点x0一个开邻域内二阶连续可微，则有二阶充分条件：若x0
% 处有gradient(f(x0))=0，gradient(gradient(f(x0)))>0成立，则x0为f局部极小点；
% 有二阶必要条件：若x0为f局部极小点，则有gradient(f(x0))=0，
% gradient(gradient(f(x0)))>=0成立。
% 称gradient(gradient(f(x)))为f在x处海瑟矩阵Hessian(f(x))，则若x0为f局部最小
% 值点，Hessian(f(x0))半正定或正定。
% 若点x处Hessian(f(x))既有正特征值又有负特征值，且gradient(f(x))=0，则称稳定
% 点x为一个鞍点。
% matlab求符号函数梯度和海瑟矩阵：
clear,clc
syms x1 x2 x3
y=x1^2+x2^2+x3^2;
x=sym([x1,x2,x3]);
h=hessian(y,x);% 求多元函数的Hessian矩阵
grady=[diff(y,x1,1),diff(y,x2,1),diff(y,x3,1)];% 求多元函数的梯度
%%  1.3.无约束不可微问题的最优性理论：
% 本小节讨论目标函数f(x)不可微的情况。
% 次梯度：设f为适当凸函数，x为f定义域中一点，向量g满足f(y)>=f(x)+g'*(y-x);，
% 则称g为函数f在点x处一个次梯度，进一步，称集合:
% SOD(f(x))={g|f(y)>=f(x)+g'*(y-x),任意y属于定义域}为f在x处次微分。
% 在本节及之后做如下定义：任意Any，属于Belongto，存在Exist。
% x为适当凸函数f的全局极小点当且仅当0 Belongto SOD(f(x))。
% 在实际问题中，目标函数不一定为凸函数，但它一般可以写成一个光滑函数与一个非
% 光滑凸函数的和，则优化问题可写为：
% min phi(x)=f(x)+h(x);其中f(x)为光滑函数(可能非凸)，h(x)为凸函数(可能非光滑)
% 复合优化问题一阶必要条件：若x0为复合优化问题局部极小点，则：
% -gradient(f(x0)) Belongto SOD(h(x0))。
%%  1.4.对偶理论：
% 本节及之后考虑更一般的约束优化问题：
% min f(x);s.t.ci(x)<=0,i Belongto I;cj(x)=0,j Belongto J
% 其中ci和cj为定义在Rn或其子集上的实值函数。
% 定义拉格朗日函数：L(x,lamda,niu)=f(x)+sum(lamdai*ci(x))+sum(niuj*cj(x))，
% 其中lamdai>=0,i=1,2,...;j=1,2,...
% 对拉格朗日函数中x取下确界可定义拉格朗日对偶函数g(lamda,niu)，即：
% g(lamda,niu)=inf L(x,lamda,niu)，则拉格朗日对偶函数g为凸函数，有原问题的
% 拉格朗日对偶问题：max g(lamda,niu);s.t.lamdai>=0,i=1,2,...
% 若原问题有最优解f(x0)，对偶问题最优解g(lamda0,niu0)，则有对偶差距：
% f(x0)-g(lamda0,niu0);若原问题约束条件线性，目标函数为凸函数，则对偶差距为0
% 若对偶差距为0，必然有：lamdai*ci(x)=0,i=1,2,...
%% 2.线性规划Linear programming
%%  2.1.matlab基础求解方法：
% 线性规划matlab求解一般形式为：
% min f(x)=f'x;s.t.A*x<=b;Aeq*x=beq;lb<=x<=ub;
% matlab基于求解器求解：
clear,clc
f=[2,3,-5];% 目标函数系数
A=[-2,5,-1;1,3,1];b=[-10;12];% 线性约束A*x<=b
Aeq=[1,1,1];beq=7;% 线性等式约束Aeq*x=beq
lb=[0,0,0];ub=[];% 自变量范围约束lb<=x<=ub
[x,fval]=linprog(f,A,b,Aeq,beq,lb,ub);% 线性规划求解器
% matlab基于问题的求解：
clear,clc
prob=optimproblem("ObjectiveSense","min");% 规定优化问题类型，求最小值为min，
% 求最大值为max
x=optimvar('x',3,'LowerBound',0,'UpperBound',[]);% 创建优化变量，添加自变量
% 范围约束lb<=x<=ub
prob.Objective=2*x(1)+3*x(2)-5*x(3);% 确定目标函数
prob.Constraints.con1=x(1)+x(2)+x(3)==7;% 确定等式约束
prob.Constraints.con2=[2*x(1)-5*x(2)+x(3)>=10;-x(1)-3*x(2)-x(3)>=-12];% 确定
% 不等式约束，不等号方向需一致
[sol,fval,flag,out]=solve(prob);% 求解优化问题
xsolve=sol.x;% 优化问题的解
%%  2.2.线性规划问题最优性检验：
% 一般线性规划问题可以转换为标准形式：
% max z=C*X;s.t.P*X=b;X>=0;一般形式中不等式约束可以通过添加松弛变量或减去一
% 个剩余变量进行转化。
% 可行基：对任意P列向量构成的满秩方阵求解，将其他变量赋值0，所有解中满足非负
% 约束的称为基可行解，基可行解对应的基称为可行基。
% 在标准形式下，求解线性规划问题流程为，求解满足非负约束的P*X=b的可行基，
% 在可行基中寻找最优解。
% 当求得当前最优解时，可进行最优性检验：
% 设将最优解变量重新排序，变成形式为x=[x1,x2,x3,...,xm,0,0,...,0];
% 求所有非基变量检验数:deltaj,j=m+1,m+2,...,n
% deltaj=cj-sum(ci*aij),i=1,2,...,m
% 若Any deltaj<0,则x为唯一最优解，若Exist deltaj>0，则可能有更优解或无最优解
% ，若Exist deltaj=0,则问题有无穷多最优解。
% matlab检验方法：
clear,clc
f=[-3,1,1,0,0];% 目标函数系数，最后两项为添加变量
A=[];b=[];% 线性约束A*x<=b，由于转化为标准型，该项为空
Aeq=[1,-2,1,1,0;4,-1,-2,0,1;-2,0,3,0,0];beq=[11;-3;1];% 线性等式约束Aeq*x=beq
lb=[0,0,0,0,0];ub=[];% 自变量范围约束x>=0，对于存在上界的变量，可以添加变量
% 转化至线性等式约束中。
[x,fval]=linprog(f,A,b,Aeq,beq,lb,ub);% 线性规划求解器
% 最优性检验
[k,~]=find(x==0);% 找到非基变量位置
delta=zeros(1,length(k));% 检验数初始化
k0=0;n=cumsum(ones(1,length(f)));n(k)=[];
for i=k'
    k0=k0+1;
    delta(k0)=f(i)-sum(f(n).*(Aeq(:,i)'));
end
%% 3.无约束优化算法Unconstrained optimization algorithm
%%  3.1.线搜索方法
% 对于无约束优化问题：min f(x)，但当前位置为xk时，一般需要确定下一步的搜索方
% 向dk和对应的步长ak。一般情况下，dk方向为当前点xk的下降方向，即：
% dk'*gradient(f(xk))<0，此时下一步的迭代点为：xk+1=xk+ak*dk;线搜索类算法关
% 键是要寻找好的方向dk和合适的步长ak。
% 本节首先讨论步长ak的选取：
% 选取ak需要满足一定的要求，这些要求被称为线搜索准则。
% Armijo准则：设dk为点xk处的下降方向，若:
% f(xk+ak*dk)<=f(xk)+c1*ak*gradient(f(xk))'*dk;，则称步长ak满足Armijo准则，
% 其中c1 Belongto (0,1)为一常数。
% 引入Armijo准则的目的是保证每一步充分下降。
% 在实际应用中可以发现当ak=0时一定满足Armijo准则，因此一般在算法设计中采用回
% 退法处理步长ak，即设定每一步迭代步长为较大的某值，在计算线搜索准则时不断减
% 小步长，直到满足准则，可以采用指数缩减:ak=gama*ak，其中gama Belongto (0,1)
% 本节介绍一种matlab结合Armijo准则和回退法的线搜索基础算法，之后的相关优化算
% 法都可以在该算法基础上进行改进。
% 线搜索基础算法Basic line search algorithm
% (该程序可单独见文件：Line_search.m)
% 求f=@(x)x(1)^2+x(2)^2;在定义域[-100,100]*[-100,100]范围内最小值
clear,clc
% 参数初始化
lb=[-1e2,-1e2];ub=[1e2,1e2];%变量范围
n=length(lb);% 变量个数
alpha0m=1/n*mean(ub-lb);% 初始步长
rand1=rand();
x0=rand1*ub+(1-rand1)*lb;% 迭代初始点
f=@(x)(x(1)-1)^2+(x(2)-2)^2;% 目标函数
% 线搜索计算过程
alpha0=alpha0m;% 初始步长
mn=300;% 迭代次数
gama2=1/mn;% 迭代次数倒数
mn1=100;% 线搜索准则判断迭代次数
x=x0;% 初始点位置
dx=(1e-6)*(ub-lb);% 求梯度时维度微元设置
gama1=0.9;% 回退法参数gama1设置
c1=1e-2;% 常数c1设置
for i=1:mn
    alpha0=(mn-i+1)*gama2*alpha0m;% 该次迭代初始步长
    tic
    alpha=alpha0;
    c_u=zeros(1,n);% 该点梯度初始化
    for bianliang=1:n% 求梯度
        c_un1=x(bianliang)-dx(bianliang);
        c_un2=x(bianliang)+dx(bianliang);
        if x(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
            % 量下限时处理
            r_un=rand();
            c_un1=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
        elseif x(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
            % 于变量下限时处理
            r_un=rand();
            c_un2=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
        end
        x1=x;x2=x;
        x1(bianliang)=c_un1;x2(bianliang)=c_un2;
        c_u(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang));
    end
    c_u=c_u/norm(c_u);% 确定该点最终梯度方向，当计算中未出现超过变量限制时为
    % 梯度方向，否则为次梯度方向。
    for j=1:mn1
        xk=x+(-c_u*alpha);% 迭代点
        for bianliang=1:n% 该步限制迭代之后的点位于自变量范围之内
            if xk(bianliang)<lb(bianliang)
                 r_un=rand();
                 xk(bianliang)=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
            elseif xk(bianliang)>ub(bianliang)
                 r_un=rand();
                 xk(bianliang)=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
            end
        end
        % 该步添加线搜索准则
        if f(xk)<=f(x)-c1*alpha*(c_u*c_u')% Armijo准则
            break
        else
            alpha=gama1*alpha;
        end
    end
    x=xk;
    toc
end
f_min=f(x);
% 在实际运算过程中，满足Armijo准则可能会使步长ak太小，此时，可以添加
% Goldstein准则：步长ak，下降方向dk满足：
% f(xk+ak*dk)<=f(xk)+c1*ak*gradient(f(xk))'*dk;
% f(xk+ak*dk)>=f(xk)+(1-c1)*ak*gradient(f(xk))'*dk;其中常数c1 Belongto 
% (0,1/2)。
% Goldstein准则可以让函数值充分下降但它可能避开了最优函数值，为此可以引入
% Wolfe准则：步长ak，下降方向dk满足：
% f(xk+ak*dk)<=f(xk)+c1*ak*gradient(f(xk))'*dk;
% gradient(f(xk+ak*dk))'>=c2*gradient(f(xk))'*dk;
% 其中c1,c2 Belongto (0,1)为给定常数且c1<c2。
% 一般情况下，以上线搜索准则均会使迭代点函数值单调下降，在实际应用中，非单调
% 算法有时会有更好的效果。
% 线搜索算法收敛性：
% 对于一般的迭代格式xk+1=xk+ak*dk，dk为搜索方向，ak为步长，当迭代过程满足
% Wolfe准则，目标函数满足下有界，连续可微且梯度-L-利普希兹连续时，当每一步下
% 降方向dk与每一步负梯度方向-gradient(f(xk))夹角thetak<pi/2时，有：
% lim(gradient(f(xk)),+inf)=0。
% 梯度-L-利普希兹连续：Any x,y Belongto Rn,Exist L>=0，有:
% ||gradient(f(x))-gradient(f(y))||<=L||x-y||
%%  3.2.线搜索方法搜索方向与步长选取分析
% 对于搜索方向dk的选取，当函数有一阶导数信息时，可以直接选择负梯度作为下降方
% 向，即：dk=-gradient(f(xk))。
% 梯度法在凸函数上的收敛性：设函数f(x)为凸的梯度-L-利普希兹连续函数，f最小值
% 存在可达，当每次迭代步长恒定且满足0<ak<1/L时，那么由梯度法
% xk+1=xk-ak*gradient(f(xk))得到的点列{xk}的函数值收敛到最优值，且在函数值意
% 义下收敛速度为o(1/k)。(本文规定：o(x)表示x的等价无穷小)。
% 当问题条件数很大，即问题比较病态时，梯度下降法收敛性会受到很大影响，此时可
% 以采用一种特殊的梯度法：Barzilar-Borwein(BB)方法：
% 引入记号sk-1=xk-xk-1;yk-1=gradient(f(xk))-gradient(f(xk-1)),求步长：
% aBB1=((sk-1)'*(yk-1))/((yk-1)'*(yk-1));
% aBB2=((sk-1)'*(sk-1))/((sk-1)'*(yk-1));
% 则可以得到BB方法两种迭代格式：
% xk+1=xk-aBB1*gradient(f(xk)); xk+1=xk-aBB2*gradient(f(xk));
% 在一般问题中，可以添加限制0<ak<amax。
% 本节给出基于非单调BB方法的matlab搜索方法：
% 非单调BB搜索算法BB-search algorithm
% (该程序可单独见文件：BB_search.m)
% 求f=@(x)x(1)^2+x(2)^2;在定义域[-100,100]*[-100,100]范围内最小值
clear,clc
% 参数初始化
lb=[-1e2,-1e2];ub=[1e2,1e2];%变量范围
n=length(lb);% 变量个数
alpha0m=1/n*mean(ub-lb);% 初始步长
rand1=rand();
x0=rand1*ub+(1-rand1)*lb;% 迭代初始点
f=@(x)(x(1)-1)^2+(x(2)-2)^2;% 目标函数
% 非单调BB搜索计算过程
alpha0=alpha0m;% 初始步长
mn=1000;% 迭代次数
e=1e-3;% 判断梯度为0点的阈值
gama2=1/mn;% 迭代次数倒数
mn1=200;% 非单调准则判断迭代次数
x=x0;% 初始点位置
dx=(1e-6)*(ub-lb);% 求梯度时维度微元设置
gama1=0.9;% 回退法参数gama1设置
c1=1e-2;% 常数c1设置
kp=5;% 非单调搜索参数设置
kpj=f(x0)*ones(1,kp);% 非单调搜索参数初始化
alpha=alpha0;% 迭代步长初始化
for i=1:mn
    alpha0=(mn-i+1)*gama2*alpha0m;% 利用BB方法计算最大步长设置
    tic
    c_u=zeros(1,n);% 该点梯度初始化
    for bianliang=1:n% 求梯度
        c_un1=x(bianliang)-dx(bianliang);
        c_un2=x(bianliang)+dx(bianliang);
        if x(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
            % 量下限时处理
            r_un=rand();
            c_un1=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
        elseif x(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
            % 于变量下限时处理
            r_un=rand();
            c_un2=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
        end
        x1=x;x2=x;
        x1(bianliang)=c_un1;x2(bianliang)=c_un2;
        c_u(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang));
    end
    G1=c_u;% 该点梯度
    if(norm(c_u)<e)% 当该点处梯度足够小时，说明找到函数驻点，跳出循环
        break
    end
    c_u=c_u/norm(c_u);% 确定该点最终梯度方向，当计算中未出现超过变量限制时为
    % 梯度方向，否则为次梯度方向。
    for j=1:mn1
        xk=x+(-c_u*alpha);% 迭代点
        for bianliang=1:n% 该步限制迭代之后的点位于自变量范围之内
            if xk(bianliang)<lb(bianliang)
                 r_un=rand();
                 xk(bianliang)=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
            elseif xk(bianliang)>ub(bianliang)
                 r_un=rand();
                 xk(bianliang)=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
            end
        end
        kpj=[f(x),kpj(2:end)];
        % 该步添加非单调搜索准则
        if f(xk)<=max(kpj)-c1*alpha*(c_u*c_u')% 非单调Armijo准则
            break
        else
            alpha=gama1*alpha;
        end
    end
    c_u1=zeros(1,n);% 迭代点梯度初始化
    for bianliang=1:n% 求迭代点梯度
        c_un1=xk(bianliang)-dx(bianliang);
        c_un2=xk(bianliang)+dx(bianliang);
        if xk(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
            % 量下限时处理
            r_un=rand();
            c_un1=r_un*lb(bianliang)+(1-r_un)*xk(bianliang);
        elseif xk(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
            % 于变量下限时处理
            r_un=rand();
            c_un2=r_un*ub(bianliang)+(1-r_un)*xk(bianliang);
        end
        x1=xk;x2=xk;
        x1(bianliang)=c_un1;x2(bianliang)=c_un2;
        c_u1(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang));
    end
    G2=c_u1;% 梯度
    % 利用BB算法计算步长
    sk=xk-x;yk=G2-G1;
    alpha=(sk*yk')/(yk*yk');% aBB1计算方法
%     alpha=(sk*sk')/(sk*yk');% aBB2计算方法
    if(alpha>alpha0)% 范围限制
        alpha=alpha0;
    elseif(alpha<0)
        alpha=0;
    end
    x=xk;
    toc
end
f_min=f(x);
% 对于正定二次函数，BB方法有R-线性收敛速度，使用BB方法的步长通常都会减少算法
% 的迭代次数。
%%  3.2.牛顿类算法分析
% 当目标函数f(x)在最优解附近二次可微时，一般可以利用牛顿类算法进行精确解求解
% 令gradient(gradient(f(xk)))=gradient2(f(xk))对于二次可微函数f(x)，在迭代点
% xk处，利用二阶泰勒展开可以得到合适的下降方向dk满足：
% gradient2(f(xk))*dk=-gradient(f(xk));，当gradient2(f(xk))非奇异时则有牛顿
% 法迭代更新格式：xk+1=xk-ak*inv(gradient2(f(xk)))*gradient(f(xk));（inv()
% 为求逆矩阵运算）。
% 经典牛顿法是收敛速度很快的算法，但它的收敛有一定条件：(1)初始点距离最优解
% 充分近；(2)海瑟矩阵gradient2(f(xk))需正定。此时我们可以通过修正牛顿法和非
% 精确牛顿法对经典牛顿法缺点进行改进，但上述牛顿类算法均需要计算函数海瑟矩阵
% 当自变量较多时计算量将会很大，实际情况下可以采用拟牛顿类算法。
% BFGS(Broyden,Fletcher,Goldfarb,Shanno)格式-Wolfe准则拟牛顿法：
% 引入记号sk-1=xk-xk-1;yk-1=gradient(f(xk))-gradient(f(xk-1))
% while 未达到停机准则 do：
% 计算方向：dk=-inv(Bk)*gradient(f(xk));
% 通过Wolfe准则找到合适步长ak>0，令xk+1=xk+ak*dk;
% 更新Bk;令k=k+1;
% 结束while
% 更新Bk方法：Bk+1=Bk+(yk*yk')/(sk'*yk)-(Bk*sk*(Bk*sk)')/(sk'*Bk*yk);
% 本节给出拟牛顿法的matlab搜索方法：
% 逆牛顿类算法Quasi Newtonian algorithm
% (该程序可单独见文件：Quasi_Newtonian_search)
% 求f=@(x)x(1)^2+x(2)^2;在定义域[-100,100]*[-100,100]范围内最小值
clear,clc
% 参数初始化
lb=[-1e2,-1e2];ub=[1e2,1e2];%变量范围
n=length(lb);% 变量个数
alpha0m=1/n*mean(ub-lb);% 初始步长
rand1=20*rand()+0.1;
f=@(x)(x(1)-1)^2+(x(2)-2)^2;% 目标函数
x_min=[1,2];% 最优解
x0=rand1*x_min;% 令迭代初始点在最优解附近
% 牛顿类算法搜索计算过程
alpha0=alpha0m;% 初始步长
mn=300;% 迭代次数
e=1e-3;% 判断梯度为0点的阈值
gama2=1/mn;% 迭代次数倒数
mn1=100;% 非单调准则判断迭代次数
x=x0;% 初始点位置
dx=(1e-6)*(ub-lb);% 求梯度时维度微元设置
gama1=0.9;% 回退法参数gama1设置
c1=1e-3;% 常数c1设置
c2=1e-3;% 常数c2设置
alpha=alpha0;% 迭代步长初始化
for i=1:mn
    alpha0=(mn-i+1)*gama2*alpha0m;% 该步迭代最大步长限制
    alpha=alpha0;
    tic
    c_u=zeros(1,n);% 该点梯度初始化
    for bianliang=1:n% 求梯度
        c_un1=x(bianliang)-dx(bianliang);
        c_un2=x(bianliang)+dx(bianliang);
        if x(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
            % 量下限时处理
            r_un=rand();
            c_un1=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
        elseif x(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
            % 于变量下限时处理
            r_un=rand();
            c_un2=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
        end
        x1=x;x2=x;
        x1(bianliang)=c_un1;x2(bianliang)=c_un2;
        c_u(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang)+eps);
    end
    G1=c_u;% 该点梯度
    if(norm(c_u)<e)% 当该点处梯度足够小时，说明找到函数驻点，跳出循环
        break
    end
    c_u=c_u/(norm(c_u)+eps);% 确定该点最终梯度方向，当计算中未出现超过变量限制时为
    % 梯度方向，否则为次梯度方向。
    for j=1:mn1
        if i==1
            dk=-c_u;
            xk=x+(dk*alpha);% 迭代点
        else
            dk=(-(B+eps)\c_u')';
            dk=dk/(norm(dk)+eps);
            if(isnan(dk(1)))
                dk=zeros(1,n);
            end
            xk=x+(dk*alpha);
        end
        for bianliang=1:n% 该步限制迭代之后的点位于自变量范围之内
            if xk(bianliang)<lb(bianliang)
                r_un=rand();
                xk(bianliang)=r_un*lb(bianliang)+(1-r_un)*x(bianliang);
            elseif xk(bianliang)>ub(bianliang)
                r_un=rand();
                xk(bianliang)=r_un*ub(bianliang)+(1-r_un)*x(bianliang);
            end
        end
        c_u1=zeros(1,n);% 该迭代点梯度初始化
        for bianliang=1:n% 求该迭代点梯度
            c_un1=xk(bianliang)-dx(bianliang);
            c_un2=xk(bianliang)+dx(bianliang);
            if xk(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
                % 量下限时处理
                r_un=rand();
                c_un1=r_un*lb(bianliang)+(1-r_un)*xk(bianliang);
            elseif xk(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
                % 于变量下限时处理
                r_un=rand();
                c_un2=r_un*ub(bianliang)+(1-r_un)*xk(bianliang);
            end
            x1=xk;x2=xk;
            x1(bianliang)=c_un1;x2(bianliang)=c_un2;
            c_u1(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang)+eps);
        end
        G2=c_u1;% 梯度
        % 该步添加Wolfe搜索准则
        if (f(xk)<=f(x)+c1*alpha*(c_u*dk')&&G2*dk'>=c2*G1*dk')% Wolfe准则
            break
        else
            alpha=gama1*alpha;
        end
    end
    c_u1=zeros(1,n);% 迭代点梯度初始化
    for bianliang=1:n% 求迭代点梯度
        c_un1=xk(bianliang)-dx(bianliang);
        c_un2=xk(bianliang)+dx(bianliang);
        if xk(bianliang)-dx(bianliang)<lb(bianliang)% 当要计算的梯度值小于变
            % 量下限时处理
            r_un=rand();
            c_un1=r_un*lb(bianliang)+(1-r_un)*xk(bianliang);
        elseif xk(bianliang)+dx(bianliang)>ub(bianliang)% 当要计算的梯度值大
            % 于变量下限时处理
            r_un=rand();
            c_un2=r_un*ub(bianliang)+(1-r_un)*xk(bianliang);
        end
        x1=xk;x2=xk;
        x1(bianliang)=c_un1;x2(bianliang)=c_un2;
        c_u1(bianliang)=(f(x2)-f(x1))/(2*dx(bianliang)+eps);
    end
    G2=c_u1;% 梯度
    % 利用BB算法计算步长
    sk=xk-x;yk=G2-G1;
    if(i==1||det(B)==0)
        B=(yk'*yk)/(sk*yk'+eps);
    else
        B=B+(yk'*yk)/(sk*yk'+eps)-((B*sk')*(B*sk')')/(sk*B*sk'+eps);
    end
    if(isnan(xk(1)))
        break
    end
    x=xk;
    toc
end
f_min=f(x);
%%  3.3.信赖域算法
% 信赖域算法和线搜索类算法都是借助泰勒展开对目标函数进行局部近似，在信赖域算
% 法中，我们直接通过求解在一个有界区域内的近似模型得到下一个迭代点，因此信赖
% 域算法在每一步迭代过程中同时选择了方向和步长。
% 对于函数f(x)，在点xk处有带拉格朗日余项的泰勒公式：
% f(xk+d)=f(xk)+gradient(f(xk))'*d+1/2*d'*gradient2(f(xk+t*d))*d;其中，t 
% Belongto (0,1)是和d有关的正数，令Bk为f(xk+t*d)海瑟矩阵近似对称矩阵，则上式
% 可变为：mk(d)=f(xk)+gradient(f(xk))'*d+1/2*d'*Bk*d;
% 此时设置deltak>0，称deltak为信赖域半径，则有信赖域
% omegak={xk+d|||d||<=deltak}，因此信赖域算法中每一步需要求解子问题：
% min mk(d),s.t.||d||<=deltak。
% 在以上子问题中仍需要确定信赖域半径deltak，我们引入如下定义来衡量mk(d)近似
% 程度好坏，以此确定deltak值：
% rhok=(f(xk)-f(xk+dk))/(mk(0)-mk(dk));若rhok接近1，说明mk(d)近似f(x)比较成
% 功，可以扩大deltak；当rhok为非常小的正数甚至是负数，此时应当缩小deltak。
% 信赖域子问题需要设计算法来快速或近似求解，可以通过迭代法、截断共轭梯度法求
% 解。由于信赖域算法，牛顿类算法都需要目标函数海瑟矩阵或其近似矩阵的信息，当
% 目标函数维度较大时会产生较大计算量，因此梯度类算法应用相对广泛。
% 本节给出基于dogleg方法的二维信赖域matlab搜索方法，该方法参考并改进了CSDN作
% 者waitingwinter原创作品，原文见：
% 原文链接：https://blog.csdn.net/waitingwinter/article/details/105305667
% 基于dogleg方法的信赖域算法Trust_region_search
% (该方法可见文件：Trust_region_search.m)
% 求解二维目标函数f = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;无约束最小值
clear,clc
% 初始化
eps=1e-8;% 梯度阈值判断，当计算点处梯度小于eps时结束算法
lb=[-100,-100];% 变量最小值
ub=[100,100];% 变量最大值
n=length(lb);% 变量个数
rand1=rand();
x0=rand1*ub+(1-rand1)*lb;% 迭代初始点
maxiter=1000;% 最大迭代次数
eta1=0.25;% rho值判断：小于eta1时说明需减小信赖域半径
eta2=0.75;% rho值判断：大于eta2时说明可以增大信赖域半径
delta=50;
delta_max=100;% 最大信赖域半径
eta=0;% 初始化eta
x=x0;% 初始化x点位置
i=0;% 初始化迭代次数
while  i < maxiter
    g = gradient(x);% 求x点梯度
    if norm(g)<= eps
        break;
    end
    h = hessian(x);% 海瑟矩阵或其近似矩阵
    p = subproblem(g,h,delta);% 利用dogleg方法求得搜索向量p
    rho =-(f(x)-f(x+p))/(g'*p+0.5*p'*h*p);% 计算rho值
    if rho < eta1 % 该部分简单判断rho值并对信赖域半径delta判断
        delta = eta1*delta;
    elseif rho>eta2 && abs(norm(p)-delta)<1e-8
        delta = min(2*delta,delta_max);
    end
    if rho>eta
        x = x+p;
    end
    i = i+1;
end
y = x;% y值为最终求得点
% Dogleg方法求解子问题
% function p = subproblem(g,h,delta)
% p_B = -inv(h)*g;% 在牛顿法搜索方向上极小值点位置向量
% p_U = -g'*g/(g'*h*g)*g;% 在梯度下降方向上极小值点位置向量
% if norm(p_B)<= delta %当牛顿法搜索方向向量小于信赖域半径时，返回该向量
%     p = p_B;
% elseif norm(p_U) >= delta % 当梯度下降方向向量大于信赖域半径时，取该方向，
%     % 步长为信赖域半径
%     p = delta/norm(p_U)*p_U;
% else % 不满足以上情况时，考虑两方向综合
%   for x = 1:10
%       disp(x)
%   end
% 
%     p_BU = p_B-p_U;
%     tau1 = (-p_U'*p_BU+sqrt((p_U'*p_BU)^2-p_BU'*p_BU*(p_U'*p_U-delta^2)))/(p_BU'*p_BU);
%     if tau1>=1 && tau1<=2
%         tau = tau1;
%     else
%         tau = (-p_U'*p_BU+sqrt((p_U'*p_BU)^2-p_BU'*p_BU*(p_U'*p_U-delta^2)))/(p_BU'*p_BU);
%     end
%     p = p_U+tau*p_BU;
% end
% end
% % 目标函数
% function f = f(x)
% f = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;
% end
% % 目标函数梯度
% function g = gradient(x)
% dg_dx1 = 2*(x(1)-1)+400*x(1)*(x(1)^2-x(2));
% dg_dx2 = 200*(x(2)-x(1)^2);
% g = [dg_dx1;dg_dx2];
% end
% % 目标函数海瑟矩阵或其近似矩阵
% function h = hessian(x)
% h_11 = 1200*x(1)^2-400*x(2)+2;
% h_12 = -400*x(1);
% h_21 = h_12;
% h_22 = 200;
% h = [h_11,h_12;h_21,h_22];
% end
%% 4.约束优化分析Constrained optimization analysis
%%  4.1.罚函数法
% 对于具有等式约束的优化问题：
% min f(x) s.t.ci(x)=0,i=1,2,...。定义二次罚函数：
% PE(x,delta)=f(x)+1/2*delta*sum((ci(x))^2);其中delta>0称为罚因子。在实际问题
% 中delta取值不能过小，否则可能出现错误情况。
% 可以编写具有等式约束的二次罚函数算法：
% 给定delta0>0，x0，k=1和罚因子增长系数rho>1，设定收敛准则，求解无约束优化问
% 题：min PE(x,delta)=f(x)+1/2*delta*sum((ci(x))^2);令：delta=rho*delta;迭
% 代求解问题。
% 对于以上算法，假设对每个delta最小值点均存在，则每次优化迭代求出的极限点列
% {xk}均为原问题的全局最优解。
% 对于具有不等式约束的优化问题：
% min f(x) s.t.ci(x)=0,i=1,2,...。定义二次罚函数：
% PI(x,delta)=f(x)+1/2*delta*sum((ci_(x))^2);其中delta>0称为罚因子，
% ci_(x)=max{ci(x),0};
% matlab函数fminsearch使用：
clear,clc
f=@(x)100*(x(2) - x(1)^2)^2 + (1 - x(1))^2;% 目标函数
x0=[-3,3];% 初始点
[x,fval]=fminsearch(f,x0);% 返回利用fminsearch找到的函数最小值点x与函数值fval
% 添加等式和不等式约束：
% 求优化问题：min f(x)=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
% s.t.x(1)+x(2)=2,x(1)-x(2)-1<=0;则：
clear,clc
delta0=2;% 罚函数因子delta初始化
rho=1.1;% 罚函数增加因子
x0=[-1,1];% 初始点
f=@(x,delta)100*(x(2)-x(1)^2)^2+(1-x(1))^2+1/2*delta*((x(1)+x(2)-2)^2+ ...
    (max([x(1)-x(2)-1,0]))^2);
delta=delta0;
sol_x=[];sol_fval=[];
for i=1:100
    tic
    f1=@(x)f(x,delta);
    [x,fval]=fminsearch(f1,x0);
    sol_x=[sol_x;x];
    sol_fval=[sol_fval;fval];
    toc
    delta=rho*delta;
end
% 从以上结果分析可知，对于fminsearch算法，其一般会找到初始点附近的极小值处；
% 当罚函数因子过小时，可能会出现结果错误。
%%  4.2.增广拉格朗日函数法
% 为了保证结果可行性，二次罚函数罚因子必须趋于正无穷，此时子问题求解可能出现
% 困难，增广拉格朗日函数法通过对二次罚函数进行修正，使得对有限的罚因子，得到
% 的逼近最优解也是可行的。对于一般约束优化问题：
% min f(x),s.t.c1i(x)=0,i=1,2,...,c2i(x)<=0,1,2,...;引入松弛变量，得到如下
% 等价形式:
%min f(x),s.t.c1i(x)=0,i=1,2,...,c2i(x)+si=0,1,2,...,si>=0;记该问题二次罚函
% 数为：p(x,s)=sum(c1i(x)^2)+sum((c2i(x)+si)^2);可以构造增广拉格朗日函数如
% 下：
% Lo(x,s,lamda,miu)=f(x)+sum(lamdai*c1i(x))+sum(miui*(c1i(x)+si))+
% delta/2*p(x,s);si>=0;其中delta>0为罚因子。
% 求解增广拉格朗日函数需要用到函数梯度信息，对于不同的目标函数需进行针对性分
% 析，对于实际问题一般更多采用二次罚函数法进行求解，本节不再进行matlab举例分
% 析。